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Logistic growth process with nonlocal interactions is considered in one dimension. Spontaneous 
breakdown of translational invariance is shown to take place at some parameter region, and the 
bifurcation regime is identified for short and long range interactions. Domain walls between regions 
of different order parameter are expressed as soliton solutions of the reduced dynamics for nearest 
neighbor interactions. The analytic results are confirmed by numerical simulations. 
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I. INTRODUCTION 

Logistic growth is one of the basic models in population 
dynamics. First introduced by Verhulst for saturated 
proliferation at a single site, it has been extended to in- 
clude spatial dynamics by Fisher [ij and by Kolmogoroff 
et. al. 0]. In its one-dimensional continuum version, 
one consider the concentration of a reactant, c{x, t), with 
time evolution that is given by the rate equation: 



dc(x, t) 
dt 



DV 2 c{x,t)+ac(x,t)-bc 2 (x,t), (1) 



where D is the diffusion constant, a is the growth rate 
and b is the saturation coefficient set by the carrying 
capacity of the medium. 

The Fisher process is a generic description of the in- 
vasion of a stable phase into an unstable region. It is 
applicable to wide range of phenomena, ranging from 
genetics (the original context of Fisher work, prolifera- 
tion of a favored mutation or gene) to population dy- 
namics, chemical reactions in unstirred reactors, hydro- 
dynamic instabilities, invasion of normal states by su- 
perconducting front, spinodal decomposition and many 
other branches of natural sciences. A comprehensive sur- 
vey may be found in recent review article by van-Saarloos 

The Fisher process ends up with a uniform saturated 
phase, in contrast with other nonlinear and reactive sys- 
tems that yields spatial structures with no underlying in- 
homogeneity. These patterns are usually related to an in- 
stability of the homogenous solution, most commonly of 
Turing or Hopf types [^|. Spontaneous symmetry break- 
ing of that type manifests itself in vegetation patterns, 
where competition of flora for common resource (water) 
induces an indirect interaction and may lead to a (Turing 
like) spatial segregation Q . 

The basic motivation of this work comes from recent 
study of non- Turing mechanism for pattern formation 
in the vegetation-water system, that yields ordered or 
glassy structures |6|. Basically, it is easy to realize that 
competition for common resource induced some indirect 
"repulsion" among agents, that may lead to spatial seg- 
regation. As an example, consider the vegetation case: 
there is a constant flow of water into the system (rain), 
and the water dynamics (evaporation, percolation, dif- 



fusion) is much faster than the dynamics of the flora. 
Now let us assume the existence of a large amount of 
flora (say, a tree) at certain spatial point. One may ex- 
pect the water density to adjust (almost instantly) to 
the tree and to equilibrate in some water profile that is 
lower around its location. The immediate neighborhood 
of the tree, though, is less favorable for a second tree 
to flourish; instead one may expect the next to grow up 
some typical distance away, reducing the water level be- 
tween them even more. This seems to be a plausible and 
generic mechanism for segregation induced by resource 
competition. These arguments may be relevant to the 
dynamics of almost any unstirred reactive system; inter- 
esting example is the process of evolutionary speciation, 
where new species may survive only " far enough" (in the 
genome space, where the spatial structure is given, say, 
by Hamming distance) from its ancestor, in order to find 
a non-overlaping biological niche. 

Surprisingly it turns out that the partial differential 
equations that describe this process (here presented in a 
nondimensionalizcd form, where w stands for water den- 
sity, b for flora and R is the "rain"): 



b(x,t) = V 2 6-^6-w6 
w(x,t) = D\7 2 w + R - w 



wb, 



(2) 



yield only a linearly stable homogenous solution. In or- 
der to get patterns one should add a cross-diffusion ef- 
fect (slowing down of the water diffusion in the presence 
of flora) that leads to Turing- like instability as in p|, 
but this is a different mechanism, and one may wander 
about the validity of the basic intuitive argument pre- 
sented above. 

Recent work 6j suggests a hint for the answer. It 
seems that a continuum and local description of a re- 
active system fails to capture the competition induced 
segregation discussed above. The continuum process is 
trying to "smear" the reactant profile, and instead of get- 
ting spatially segregated structure of large biomass units 
(" trees" ) it favors homogenous profile of " grass" cover- 
ing all the area. In 0] a biomass unit was allowed for 
long time survival only if it exceeds some prc-dctcrmincd 
threshold, and simulation of the system reveal an imme- 
diate appearance of spontaneous segregation and stable 
patterns. 

Similar situation appears, presumably, in the process 
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of bacterial colony growth where the food supply is lim- 
ited. As noted by Ben-Jacob et. al. 0, spatial seg- 
regation and branching are induced by the competition 
of bacteria for diffusive food. A communicating walkers 
model is used by these authors to simulate the branching 
on a petri-dish, where continuum equation dictates the 
food dynamics while the individual bacteria are discrete 
objects. The discreteness of bacteria adds some weak 
threshold to the system and induces segregation. Note, 
however, that the model admits weak dependence of the 
diffusion on the local bacterial density at the boundaries 
of the colony. 

In this work I consider the one-species analogy of the 
competition problem, namely, a logistic growth with non- 
local interactions, where the carrying capacity at a site 
is reduced due to the presence of "life" in another site. 
Nonlocal competition has been recently considered by 
Sayama et al. @ and by Fuentes et al. [£| . Both groups 
uncover the possibility of spontaneous symmetry break- 
ing and patterns, depending on the the strength and the 
smoothness of the "weight function" that controls the 
nonlocality. It seems that nonlocal interactions are not 
simply an effective model obtained by integrating out 
the fast degrees of freedom; rather, it incorporates some 
nonlinear effects (like the threshold) and allow for lin- 
ear instability that manifests the intuitive "competition 
induced segregation" argument. 

Sayama et al. Q deal with a two dimensional model 
of population dynamics, with no diffusion term. Both 
the local growth term and the carrying capacity at a site 
depend (not in the same way) on the population of neigh- 
boring sites; in a crowded neighborhood the growth term 
becomes larger (due to offspring migration) while the car- 
rying capacity decreases as a result of long range competi- 
tion. The conditions for an instability of the homogenous 
solution have been found analytically and demonstrated 
numerically for a "stepwise" weight function (taking as 
the effective neighborhood the average density inside a 
prescribed radius around the site). It was also pointed 
out that a gaussian weight function yields no instability. 

In the numerical work of , a one dimensional realiza- 
tion of diffusing reactants has been considered, equivalent 
to Fisher equation with non-local interactions. Again it 
was shown that a stepwise weight function may lead to in- 
stability while a gaussian weights lead to stable homoge- 
nous solution; the authors proceed to consider interme- 
diate weight functions that interpolate between gaussian 
and a step function. 

As in any case of spontaneous symmetry breaking, the 
system falls locally into one of the "minima" of the or- 
der parameter, and typically domains are formed. These 
domain walls determine the low lying excitation spec- 
trum of the system, as their movement is "smooth": if 
the broken symmetry is continuous the resulting Gold- 
stone modes may destroy the long range order at finite 
temperature, and the same is true for the domain walls 
if discrete symmetry is broken. Although we are deal- 
ing with an out of equilibrium system, one may guess 



that the response to small noise is determined by these 
domain walls. 

The goals of this work are twofold: in the next section 
I will try to give more comprehensive discussion of the 
instability condition, with and without diffusion, and its 
dependence on the weight functions: it turns out that it 
depends on the minimal value of its Fourier transform. 
The third section is devoted to the appearance of topolog- 
ical defects in the segregated phase. Finally in section IV 
some discussion and possible implications are presented. 

II. INSTABILITY CONDITIONS 

The model is a one dimensional realization of long-range 
competition system on a lattice (with lattice spacing Iq) 
and the continuum limit is trivailly attained at Iq — > 0. 

In the generic case of diffusion and non-locality the 
time evolution of the reactant density at the n-th site, 
c„, is given by: 



Q t = 2c„(t) + c„ + i(t) + c„_i(t)] + ac n (t) 

OO 

- b& n {t)-C n (t)J2lr[Cn+r(t) + C n - r (t)]. (3) 
r=l 

where D is the diffusion constant and a, 6, 7 are the corre- 
sponding reaction rates. The definition of dimensionlcss 
quantities, 

t = at, c — bc/a, r y r = 7 y r /b, D = —pj. (4) 

(the new "diffusion constant" is D = W 2 /Iq, where 
W = ^jD/a is the width of the Fisher front) provides 
the dimensionless equation, 

dc n 

— — = D[-2c„ + Cn+i + C„_i] 

OT 

+ C n f 1 - Cg - -y r [c n +r + Cn-r]^ , (5) 

that may be expressed in Fourier space [with = 
En c n e lknLo } as, 

A k = a k A k -^2/3 k ^ q A q A k _ q , (6) 

Q 

where 

a k = 1 - 2D[1 - cos(kl )} (7) 

OO 

l\ = l + 2Y,lrCos(rkl ). (8) 

r=l 

As c„ is positive semi-definite, A$ is always "macro- 
scopic". Any mode is suppressed by Aq, and for small 




FIG. 1: Maximal amplitude of c n differences [max(c n ) — 
min(c n )] (circles) for a sample of 1024 sites (periodic bound- 
ary conditions) and the predicted difference according to Eq. 
Illl l for nearest neighbor interaction with 7 = 0.8. (dashed 
line). The agreement is up to the numerical error. 



7 r one expects only the zero mode to survives 10]. If, on 
the other hand, j r increased above some threshold, bi- 
furcation may occur with the activization of some other 
k mode(s), and the homogenous solution becomes un- 
stable. This is the situation where patterns appear and 
translational symmetry breaks. 

To get a basic insight into the problem, let us consider 
the case with no diffusion (£) = 0, a* = 1). Eq.JHJ) 
becomes, 



(9) 



and division by c„ yields, for the steady state, the linear 
equation Q.c — y, where Q is a circular matrix, c is the 
vector of c„-s and y = (...1, 1, 1, 1, The sum of the 
elements of any row of Q is the same, so the homogenous 
state (scalar multiplication of y) should be an eigenvec- 
tor. On the other hand, if Q is non-singular it must admit 
a full set of mutually orthogonal eigenvectors. Only the 
constant eigenvector of Q has nonvanishing projection on 
y, so the only positive definite, non diverging steady state 
(c„ = c n > Vn.) solution is the homogenous one, 

Cn = 1/Pq. 

As implied by Eq. © , the homogenous steady state is 
unstable iff, for some k, 0k < 0. In that case bifurcation 
occurs, and the new steady state is a combination of the 
zero mode and the k mode with equal weights Aq = Ak = 
1/(00 + 0k). 

The function 0k, k € [0, ft /la], is discrete for finite 
systems and becomes continuous at the thermodynamic 
limit. If /3k never crosses zero there is no bifurcation and 
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FIG. 2: A typical "soliton" of length L — 20, an outcome of 
forward Euler integration of Eq. (I13P on 1024 lattice points 
with periodic boundary conditions and random initial condi- 
tions at 7 = 0.505 (just above the bifurcation). There is a 
perfect agreement with the theoretical prediction, Eq. 1211 . 
up to the accuracy of the numerics (here, 4-5 significant dig- 
its) 



the homogenous solution c„ = l/0o is stable. The results 
for few types of interaction ranges, with the critical value 
7i (where the instability occurs) , and k c (the first excited 
mode), are summarized in Table I. 

It is interesting to note that these expressions may be 
generalized to yield a full, period doubling type, instabil- 
ity cascade. The m-th instability involves 2 m modes, and 
the steady state is 1/ /3k, where the sum runs over all 
the "active" modes. The condition for the m + 1 bifurca- 
tion [activation of another 2 m+1 modes] is the existence 
of a wavenumber q such that 0q-k < 0, with the 
sum runs, again, over all 2 m active k-s. There are, how- 
ever, some obstacles for the implication of these solutions 
above the first bifurcation. Degeneracies in (3k (e.g., for 
7r = <5r,4, both Hq = 7r/4 and Mq = 37r/4 are minima), 
and solitons between different stable phases (described 
below) may blur the native state. In this letter, though, 
(3k is used only for the first, pattern-forming, instabil- 
ity criteria, and the details of the emerged structure are 
presented just for nearest-neighbor interaction. 

Once diffusion is added to the system, its features 
changes, but not so much. The homogenous state is still 
characterized by c„ = l//?o and the first pattern forma- 
tion instability appears when some k mode satisfy: 



0k < -20 O D[1 - cos(kl ) 



(10) 



Above this instability, the amplitudes of the modes are 
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type 


7r 


Pk 


Instability condition 


exponential 
quadratic 

step 
gaussian 


(7i) r/5 

r 1 ' — V 

r > p 

(7i) (rVCT2) 


sinh(\ln{-yi)\/£) 
cosh( Zra(7i) cos (fcio) 

l + 2 7l [^ T«o + («o) a ] 

-„ / ka /<-)„ i i \ \ 

if^»io,~!/^"e, P ( 4 ^) 


no instability 

7i = Jy, fc c = tt/^o 
at large p k c = 

no instability 



TABLE I: The function /?;, and the instability conditions for various types of nonlocal interactions. The results for the Gaussian 
case are in the continuum approximation. 



not equal, 



A = 



(3 + (3 k 



A k = 



MPo+Pk) 2 



(11) 



and there are no zeroes of c n . This result fits perfectly 
with the numerical data presented in figure Q). Again 
the m-th instability involves the activation of 2 m modes, 
Although the stability analysis is more complicated. 

The question of pattern instability is thus translated to 
the determination of the minimal value of the Fourier co- 
efficient of the weight function (or the "weight series" 7 r ). 
If the minimal value is smaller than some prescribed num- 
ber (zero if there is no diffusion) instability takes place 
and patterns emerge. Unfortunately I am not familiar 
with a general theorem that sets bounds on the minimal 
value of the Fourier coefficient of a function based on its 
"smoothness", or other analytic properties, so any case 
should be considered separately, with the generic exam- 
ples given in Table I. 



III. DOMAIN WALL STRUCTURE 

Above the pattern formation threshold generic initial 
conditions fail to yield perfect " lattice" , as different do- 
mains reach saturation with different "phases". These 
domains are connected by soliton-like solutions of the 
time independent equation in the following sense: any 
stable solution (c(x, t) — 0) should satisfy (in the contin- 
uum limit ), 



D 



d 2 c{x) 
dx 2 



= c{x) - c(x) / f(x - y)c(y)dy, (12) 



thus it looks like a trajectory of mechanical particle (with 
mass D) in a nonlocal potential, with x as the "time". A 
domain wall is a finite size structure, so it must connect 
fixed points of this fictitious dynamics, i.e., a domain wall 
corresponds to heteroclinic orbit. In this section we con- 
sider these " solitons" and look for their shape and size at 
different conditions. In order to simplify the discussion, 
only the nearest neighbor case is considered, both with 
and without diffusion. 



With no diffusion and n.n. competition, © takes the 
form: 



den 
dr 



Cn [1 - C„ - 7<Ai+l + C n _l)] 



(13) 



The uniform solution, in this case, is c — - — r- 
' ' 1+27 • 

1 for odd n and c 



' and the 




nonuniform solution is, either c 
for even, or vice versa. Stability analysis shows that the 
uniform solution becomes unstable at 7 C = 1/2, and the 
zero-one phase is stable above this value. One may ex- 
pect, though, to see a jump from homogenous to pat- 
terned (zero-one) phase at j c . However, if the initial 
conditions are taken from random distribution, there is 
a chance for a domain wall between two regions, as indi- 
cated by the numerical results presented in Figure J2J- 

Clearly, such a soliton should be a solution of the 
"map" 



1 



C n +1 



7 



(14) 



of course, 01010101... (odd 0-s) and 101010101... (even 
0-s) are already solutions of this equation. We are look- 
ing for the solution that connect these two fixed points. 
Such a trajectory begins in, say, 010101010 state, but 
then after the zero it gives not 1 but x\. The dynam- 
ics now continue in a different trajectory, but x\ should 
be selected such that after L steps of the map (for do- 
main wall of size L) the 101010 solution is rendered. In 
a matrix form, the condition that determined xf {x\ for 
a given L) is: 



,L I - 



7 

1 





7 


1 



- L 




















= M L 













(15) 



Where we assume symmetry of the soliton, so L must be 
even. In other words, the condition that determine x\ 
is: 
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" _l 
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-1 
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(16) 
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FIG. 3: Solitons for nn interaction with 7 = 0.505 and D = 
0.001 (circles) D = 0.0012 (heavy line) and D = 0.00124 
(line). D c = 0.00124378 



Diagonalization of M. is given by the matrix S: 



where 





S^MS = 


V 


(17) 




1 







V = 


-e" je 





(18) 












and 8 = arctan(^/47 2 — 1) = arccos(^). 

The eigenvalue problem 1)16(1 may be written in terms 
of the diagonal matrix, 



L c c r 13 CQS(— + tp) 

X\ — J11 — ^12 TTfl r 

T23 cos(4f +rj) 



(19) 



and using 6*12 = Sl z = r^e^ and S22 — ^23 = r 23e l?? it 
implies that, 



1 



2 7 + l 




i-r £ _!(i/2 7 ) 
i-r L+1 (i/2 7 ) 



(20) 



where T n {x) is the n-th Chebyshev polynomial of the 
first kind. After determining x\ the same method may 
be used to derive a general expression for all the elements 
of the size L soliton (x^, where 1 < m < L): 



27+I 



i + (-ir 



/ 2 7 [l-T 2m _ 1 (l/27)] 



27-I 



-(-l) m ^2^Vl + TWl/27) 



(21) 



and it fits perfectly the numerical experiment presented 
in Fig. see captions). 

The above analysis gives the shape of a soliton for any 
prescribed length L,but simulation indicates that only 
one soliton size L is selected for any set of parameters, 
and its length diverges as 7 approaches its critical value. 
Looking carefully at the solutions (|21|) one realizes that 
all other possible solitons admit values for some of the 
x m -s that are either negative or larger than one, so these 
options are unphysical (negative) or unstable to small 
perturbations. 

The actual £(7) may be forecasted by a rough ar- 
gument based on a continuum approach. Defining the 
local deviation from the one-zero solution, c„ = 1 and 

Cn±l = 0, 



Cn±l — §n±l C n — \ — S n 

and plugging it into Eq. ((Tljl gives, 



Sn+1 + Sn-1 



7 



(22) 



(23) 



Subtracting of 2S n from both sides and taking the con- 
tinuum limit (i.e., assuming that the changes in S from 
site to site are small compared to Iq, here taken to be 
unity) one gets, 



Ac 



1 + 2e 



S(x) 



(24) 



with e = 7 — 7 C goes to zero at the transition, so 1 + 2e « 
1. The solution of Eq. (|24ll that satisfies the boundary 
conditions 6(0) — together with S(L) = 1 is 



5(x) 



sin(2^/e)L 



(25) 



This expression fails to converge smoothly to the "back- 
ground" ordered 010101 phase (at x=0 it has a finite 
slope), and is also asymmetric. Put it another way, there 
are no nontrivial heteroclinic orbits for a parabolic po- 
tential. On the other hand, close to the transition, where 
the size of the domain wall is large, it seems that it should 
fit a solution of the continuum approximation. The only 
way out is to pick a soliton size L such that the continuum 
equation admit no solution at all, i.e., 



L = 



17 



(26) 



Such a choice forces us back to the discrete equation ltT4"|) 
and its solution . This argument turns out to give the 
right length of the stable soliton in the large L (e — > 0) 
limit, as shown in Fig. J3J|. 

Let us consider now the domain walls for the finite 
diffusion case. As seen in Fig. (0, there are also soli- 
tons for the finite diffusion case, but now they admit 
tails that asymptotically looks like 5 ~ exp(—x/t;), and 



£ diverges at the transition [e.g., at D c = 



27-1 
4(27+1) 



for 



6 




FIG. 4: Soliton size (L) as a function of (e = 7 — 

7c), for a nearest-neighbor interaction without diffusion. The 
circles are the results of a numerical simulation and the line is 
n/(2y/e). In the inset, the characteristic length of the soliton 
tail, £, is plotted against 1/y/rj for 7 = 0.505 for the solitons 
shown in Fig. ©. The circles come from the numerics and 
the line is the best linear fit that give a slope of 0.198, to be 
compared with 1/V32 = 0.177. 



nearest neighbor interaction]. Defining a vectorial "or- 
der parameter" according to the larger c values, soliton 
solution interpolates between (1,0) (larger c on the odd 
sites) to (0, 1) (even sites), and its shape is given by the 
saddle point solution of the appropriate dynamics. Al- 
though the determination of its full shape is difficult, it 
is possible to determine £ by linearizing around one fixed 
point. For small deviations, c„ = Aq + A n n — S n and 
Cn±l = A — A^ji Q + 5 n ±i and Eq. (|T4"|l yields the two 
coupled linear equations for odd and even n-s: 



(2s + 2d 1 -l + 2D)5^ n + (D- 7S )(«S°f 1 + 6%*) = (B7) 
{2d + 2s 7 - 1 + 2D)8° n M + {D- 7 d)(C+T + CT) = 0. 



where s = Aq + A n n and d = Aq — A^n^. These cou- 
pled equation may be solved with the ansatz 5 even ~ 
ai exp(— x/£), S odd ~ a2 exp(— x/£) to give, 



, / 1 / (4 + 64L> 2 - 32Z))7 2 + (2QD - 8D 2 - 4)7 + 1 - 20D 2 

arccosh —4/ ; ; . 

V 2 V D( 7 - D - 2jD) 



, -1 



(28) 



As the diffusion constant approaches its critical value, 
D = D c — 77, £ diverges like 1/^32^. This prediction 
is tested in the caption of Fig. Q against the numerics 
and there is reasonable quantitative agreement, given the 
difficulties in getting reliable numerical accuracy for the 
slope of the logarithmic tail of the soliton. 



IV. CONCLUDING REMARKS 

The model of logistic growth with long range interac- 
tion term may serve as a generic, minimal model for com- 
petition for common resource and pattern formation in 
excited media. In this paper this model has been analyti- 
cally discussed, with two main outcomes. First, a general 
scheme for the identification of the pattern forming in- 
stability has been presented, along with explicit results 
for few common cases. Second, the defected solutions for 
random initial conditions has been presented and ana- 
lyzed, and their characteristic length that diverges at the 
transition is calculated. 

The patterned solutions (like the 010101 for n.n. inter- 
actions) are stable against small perturbation (like a low 



amplitude white noise). If instead of ...01010101.. one 
have ...01010(0.9)01.. the 0.9 site relaxes to 1. The do- 
main walls, on the other hand, are much less stable, and 
their density and dynamics has to be strongly effected by 
an external noise. This is very much like the situation in 
magnetic system, where the response functions of the ma- 
terial are basically determined by the domain walls and 
not by the "bulk". In magnetic systems, however, one 
may define the state of the system at finite temperature 
ad a minimum of the free energy function and consider 
the effect of noise simply as temperature increase. The 
situation seems to be different for the long range compe- 
tition model: no simple Liapunov function exist for this 
system, and the steady state is not derived from some 
variational principle. In spite of that it is plausible to as- 
sume that the defected solutions determine the response 
function of the segregated phase, and maybe an effective 
dynamical equations for the solitons may be set up to 
give an approximate Liapunov function for this system. 
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